Crystal Nucleation in a Supercooled Liquid with Glassy Dynamics 
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In simulations of supercooled, high-density liquid silica we study a range of temperature T in 
which we find both crystal nucleation, as well as the characteristic dynamics of a glass forming 
liquid, including a breakdown of the Stokes-Einstein relation. We find that the liquid cannot be 
observed below a homogeneous nucleation limit (HNL) at which the liquid crystallizes faster than 
it can equilibrate. We show that the HNL would occur at lower T, and perhaps not at all, if the 
Stokes-Einstein relation were obeyed, and hence that glassy dynamics plays a central role in setting 
a crystallization limit on the liquid state in this case. We also explore the relation of the HNL to 
the Kauzmann temperature, and test for spinodal-like effects near the HNL. 

PACS numbers: 64.60.My, 64.60.Q-, 64.60.qe, 61.20.Lc 



A long-standing question regarding the liquid state 
concerns its ultimate fate when supercooled below the 
equilibrium freezing temperature [![. In his seminal 1948 
work on the glass transition, Kauzmann pointed out that 
many supercooled liquids are headed for an "entropy 
catastrophe" as the temperature T decreases |2|. That 
is, the liquid entropy decreases so rapidly with T that 
unless crystallization or the glass transition intervenes, 
the entropy will become negative below a finite temper- 
ature Tx . Since 1948, two major scenarios have emerged 
for avoiding the entropy catastrophe based on our un- 
derstanding of glasses: that the liquid terminates at Tk 
in an "ideal glass" state; and that a "fragile-to-strong" 
crossover occurs in which the entropy changes its In- 
dependence so as to remain non-zero for all T > [H, Q • 

Yet it is often crystallization, rather than the glass 
transition, that terminates the liquid on cooling prior 
to Tk- In experiments, a homogeneous nucleation limit 
(HNL) is often encountered below which crystallization 
is, in practice, unavoidable [l[. In theoretical and simu- 
lation studies, a kinetically-defined HNL has been iden- 
tified as the T at which the time for crystal nucleation 
T n becomes comparable to the structural relaxation time, 
T a @, H, 0, [f|. Below this HNL, the supercooled liquid 
ceases to be observable because it nucleates before it can 
equilibrate. Other recent simulation studies report the 
occurrence of a thermodynamically-defined spinodal-like 
limit for crystal nucleation, at which the free energy bar- 
rier to nucleation AG* decreases to zero[(| 0, [||. 

Kauzmann himself proposed that crystallization be- 
comes inevitable in supercooled liquids prior to Tk 0, 
and recently Tanaka has reexamined this scenario, 
through an analysis of both classical nucleation theory 
(CNT) and the behavior of glass forming liquids [3(. 
While CNT forms the basis of much of our understand- 
ing of nucleation, by itself it seems to predict no limit 
on supercooling [1| . In CNT the nucleation time is given 
by t„ = K exp(AG* / RT) , where R is the gas constant, 
and the kinetic prefactor K contains a factor of D -1 , the 



inverse of the liquid diffusion coefficient. Assuming that 
D^ 1 is proportional to t q , a system obeying CNT will 
always satisfy t„ > t q , i.e. the equilibrium liquid will be 
observable prior to nucleation at all T. 

In this context, Tanaka pointed out that a liquid can 
obey CNT and exhibit a HNL [9j if there is a breakdown 
of the Stokes-Einstein (SE) relation The SE relation 
asserts that Drj/T ', where rj is the viscosity, is indepen- 
dent of T. Here we use t q as a proxy for 77, since both 
quantify collective structural relaxation. The violation 
of the SE relation is a ubiquitous feature of supercooled 
liquids, in which Dr a /T is found to grow rapidly as T 
decreases, demonstrating that the relaxation time associ- 
ated with diffusion (a local process) increases faster than 
r a (a global process). Since in CNT t„ is controlled by 
D, and not by r Q , SE breakdown makes it possible for 
T n and t q to become comparable at sufficiently low T, 
inducing a HNL. 

Tanaka's analysis is significant because it predicts that 
the physics of a glass forming liquid (SE breakdown) is 
crucial to the origins of the phenomenon (unavoidable 
crystallization at the HNL) by which the entropy catas- 
trophe is avoided. Tanaka showed that experimental data 
on a metallic liquid alloy is consistent with his interpre- 
tation Q. Cavagna and coworkers have come to similar 
conclusions by incorporatingSE breakdown into CNT 
through viscoelastic effects [!|. In this Letter we use 
computer simulations to identify a HNL in a deeply su- 
percooled liquid, and show that the nature and location 
of this limit is indeed strongly influenced by the presence 
of glassy dynamics. 

Our results are based on molecular dynamics simula- 
tions of the BKS model of silica (Si02) [10| . Our system 
consists of 444 Si atoms and 888 O atoms at fixed volume, 
with simulation parameters the same as in Refs. ITlL [T^] . 
This model liquid is a well-studied glass- former [13j , and 
yet also crystallizes to stishovite on time scales accessible 
to simulation for T < 3200 K at density 4.38 g/cm 3 [TO]- 
Under these conditions the liquid is deeply supercooled; 
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FIG. 1: (a) Behavior of r a , T n and t^ b as a function of T. The 
model functions r,^ FT (solid line), t^ ft /R (dashed line) and 
£ T a FT (dot-dashed line) are also shown, (b) T-dependence 
of r n /r a and r n /r„ E . The model functions £ (solid line) and 
£R (dashed line) are also shown. 



at this liquid density the liquid-stishovite coexistence 
temperature is greater than 6000 K [ill ]. 

We quantify the dynamical properties of the liquid by 
evaluating D and r a at density 4.38 g/cm 3 for several 
T from 5000 to 2900 K. D is evaluated from the mean 
square displacement, while r Q is defined as the time con- 
stant in a fit of a stretched exponential exp[— (i/ T a) ] 
to the decay of the intermediate scattering function at a 
wavenumber corresponding to the first peak of the struc- 
ture factor. Both quantities are computed for the Si 
atoms only. For each T > 3200 K, these properties are 
evaluated from a microcanical simulation run, starting 
from a well-equilibrated initial configuration. 

For T < 3200 K, crystallization occurs spontaneously 
on our computational time scale. In this regime we seek 
to evaluate D and r a for the liquid as well as to quan- 
tify the crystal nucleation kinetics. To this end, for each 
T < 3200 K, we conduct 200 runs initiated from distinct 
configurations equilibrated at 5000 K, and quench the 
system by applying the Berendsen thermostat [l4[ (with 
a time constant of 1 ps) to decrease T to the target value. 
Each run continues until it crystallizes, as detected by a 
significant drop in the potential energy. The nucleation 
time for a run is taken as the latest time at which the 
system contains no crystalline particles [l5| • These times 
are averaged over the 200 runs to give the mean nucle- 
ation time T n ; uncertainties are computed as the standard 
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FIG. 2: (a) VFT and (b) AG plots of D/T (circles) and l/r Q 
(squares). Filled symbols are data generated for T > 3200 K 
from single simulation runs. Open symbols are data found 
using the 20 longest runs from the ensemble of 200 conducted 
at each T < 3200 K. Fits of a straight line to each data set 
are also shown. In (a) T = 2254 K. In (b) the data for S c 
from Fig. |3ja) is used to evaluate 1/TS C . 



deviation of the mean. 

To estimate D and r a for each T < 3200 K, we iden- 
tify those 20 of the 200 runs that remain longest in the 
liquid state before crystallizing. At three T (3100, 3000 
and 2900 K) we find that the transient behavior asso- 
ciated with the quenching procedure is completely re- 
moved by discarding the initial 2 ns of each run. The 
liquid state properties (e.g. energy, pressure) are sta- 
tionary thereafter, and we find D and r a for each run 
from this stationary time series. Averaging over the 20 
runs, we compute the mean value of D and r Q , and their 
respective uncertainties as the standard deviation of the 
mean. In this work, we restrict our attention to those 
T at which we can confirm that liquid equilibrium is es- 
tablished on a time scale much less than the nucleation 
time (i.e. r a << r n ) in order to ensure that nucleation 
events do not interfere with the accurate evaluation of 
liquid properties. 

It is important to note that t„ reported here is the 
system nucleation time, which depends on the system 
volume V as T n = (JV)~ 1 where J is the nucleation rate 
per unit volume 1]. Hence r„ for larger systems than 
ours will be smaller by a factor of V /V, where V is 
our system volume. In simulations, as in experiments, 
small systems are often exploited in order to increase the 
nucleation time and thus allow examination of deeply 
supercooled liquid states. As shown below, our system 
size is small enough to allow us to reach a range of T in 
which both nucleation and glassy dynamics occur; and 
yet large enough to easily accommodate crystal nuclei of 
critical size. 

In Fig. QJa) the T-dependence of r„ is compared to 
that of T a . Previous studies indicate that liquid equi- 
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FIG. 3: (a) S c as a function of T (thick solid line). The dashed 
line is an extrapolation based on the finding in Fig.[2]that both 
the VFT and AG relations are satisfied. That is, we solve 
Aexp[B/(T-T )} = C exp[E / (T S c )} for S c , where the fitting 
parameters A, B, T , C and E are taken from the fits to D/T 
in Fig. [2] This curve models the case where S c approaches 
an ideal glass transition at Tk = T . (b) T dependence of 
the SE ratio Dr a /T, normalized by c D , the value of Dr a /T 
at 5000 K. The solid line is the model function R. 



librium is established on a time scale of between 10 and 
20 times r Q [H, [T^- At T = 3100, 3000 and 2900 K, 
we find that r n remains greater than an order of magni- 
tude larger than r Q , indicating that liquid equilibrium is 
well established prior to crystal nucleation. However, the 
gap between r„ and r a is closing rapidly as T decreases. 
Fig. QJb) shows that the ratio T n /r a is decreasing with 
T in a manner suggesting that this liquid system reaches 
a HNL in the vicinity of 2800 K, where nucleation will 
on average occur faster than equilibrium measurements 
of liquid properties can be made. Consistent with this 
behavior, we have attempted simulations at 2800 K, but 
find that a significant fraction of the 200 runs nucleate 
during the initial transient associated with the quench- 
ing procedure, complicating the evaluation of r„, as well 
as of D and r a . Empirically, we also note that in the 
T range studied here r n jr a fits well to an exponential 
function (referred to below as £), as shown in Fig. [T^b). 

On approach to the HNL, we find that the liquid be- 
haves as a fragile glass former, in that D/T and r a are 
modeled well by the Vogel-Fulcher-Tammann (VFT) ex- 
pression Aexp[B/(T - T )}, where A, B and T a are 
fitting parameters [Fig.^a)]. The dynamical divergence 
temperature T a serves as an estimate of Tk- We obtain 
T = 2270 and 2237 K for D/T and r Q respectively, and 
in Fig.[5]Ja) we set T a to the mean of these values, 2254 K. 
The VFT fit to r Q , denoted r^ FT , is shown in Fig. [Ha). 

The fragile nature of the liquid is also demonstrated 
in the T-dependence of the configurational entropy S c 



FIG. 4: AG(n) /RT as function of n for T = 2900, 3000, 3100, 
3200 and 3300 K (data sets from bottom to top). Solid curves 
are fits to the CNT prediction AG(n)/RT = -an + bn 2/3 , 
where a and b are fitting parameters. Inset: AG* /RT and n* 
as a function of T. 



[Fig-Et a )]- S c quantifies the entropy associated with the 
number of distinct basins of the potential energy land- 
scape explored by the liquid at a given T. We evaluate 
S c via an analysis of the inherent structure energy of the 
liquid, as described in detail in Ref. [12}. We find that 
S c decreases rapidly as T decreases, characteristic of a 
fragile glass former. Using S c we also find, in common 
with many glass forming liquids, that D/T and r Q satisfy 
the Adam-Gibbs (AG) expression Cexp[E/(TS e )}, 
where C and E are fitting parameters [Fig. [2jb)]. The 
fact that D/T and r a conform to both the VFT and AG 
relations at all T in Fig. [2 is a validation of our method 
for finding the equilibrium liquid behavior in the low T 
range, where nucleation also occurs. 

A further signature of glassy dynamics is shown in 
Fig. EJb), which gives the SE ratio Dr a /T normalized 
by its high-T value c a . As T decreases, the SE relation 
breaks down, and at the lowest T the characteristic time 
scale for structural relaxation is nearly 10 times larger 
than that associated with the diffusion process, compared 
to high T. We use the VFT fits to D/T and t q to eval- 
uate a model function (denoted R) for DT a /Tc , shown 
in Fig.E^b). 

Next we seek to quantify the role played by SE break- 
down in the occurrence of the HNL. To do so, we estimate 
what r a would be if the SE relation were obeyed. That 
is, we define t^ e as the value of r a computed using the 
data for D via the SE relation, i.e. rf E = c T/D. If we 
assume that r n depends on D, but not on r Q , as predicted 
by CNT, then r„ can be compared to rJ E to test for the 
occurrence of a HNL if the SE relation were obeyed. 

We find that t~ e is as much as an order of magnitude 
smaller than r a , and that r n and rJ E are more widely 
separated, and converging less rapidly, than t„ and r a 
[Fig. Ufa)]. In Fig.QJb) we plot the ratio r„/r^ E . Noting 
that t„/t^ e = (t„ /Tq.) (DT a / Tco) , we model r„/rf E as 
the product of the fitting functions £ and R. In Fig.flja), 
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we model rf E as rJ FT /i?, and r„ as £t^ ft , to show 
the behavior consistent with the modelling presented in 

Fig.njb). 

Comparison of r„/r Q and t„/t^ e in Fig.rjTb) suggests 
that if the SE relation were obeyed, the HNL would 
be shifted to lower T by at least several hundred de- 
grees K. Since r„/r Q and t„/t^ e differ by a factor of 
Dr a /Tc , which appears to diverge as T — > T Q , the gap 
between r„/r a and t„/t^ e will grow rapidly as T de- 
creases. So long as T n /r a decreases with T more slowly 
than DT a /Tc Q increases, it is possible for the liquid to 
remain observable (i.e. satisfy r a << t„). For example, 
if the exponential fit to T n /T a holds for T < 2900 K, our 
data allow the possibility that a liquid obeying the SE 
relation remains observable approaching T . 

In sum the above results illustrate a case in which SE 
breakdown plays a key role in setting a supercooling limit 
on the liquid state due to crystallization. This interplay 
of glassy phenomena and crystal formation is a realiza- 
tion of Kauzmann's original proposal for avoiding the en- 
tropy catastrophe at Tjf, and is entirely consistent with 
Tanaka's recent analysis Q. In our system, Tjc cannot be 
approached because of a HNL at approximately 2800 K. 
Yet, in the absence of glassy dynamics (in the form of 
SE breakdown), lower T would be accessible, including 
perhaps the region near Tk- 

We note that the effect of system size on r„ explained 
above will not qualitatively change our conclusions. De- 
creasing V, for example by a factor of 10 (from 444 to 44 
molecules), will increase r„ by a factor of 10 and lower the 
HNL by approximately 200 K. This would place the HNL 
at approximately 2600 K, still well above T a — 2254 K. 
It would be difficult to justify any further reduction in 
V that would not introduce significant, unphysical finite- 
size effects. The HNL studied here is thus located near 



the lowest possible T at which such a limit can be ob- 
served in this system, and yet is always well above T a for 
all V. 

Finally, we also examine the thermodynamic aspects of 
the nucleation process, to test if the kinetically-dcfincd 
HNL is related to a spinodal-like thermodynamic limit. 
We show in Fig. [5] AG(n), the work of formation of crys- 
talline clusters of size n for several T. Our procedure for 
computing AG(n) is the same as that used in Ref. [TTI j . 
The number of molecules in the critical nucleus n*, as 
well as AG* are both decreasing as a function of T (in- 
set, Fig. QJ. At the same time, both quantities remain 
finite in the range of T studied here, and the shape of 
AG(n) remains consistent with CNT. Further, in nucle- 
ation influenced by a spinodal, it is expected that the 
critical nucleus becomes ramified and/or that the nucle- 
ation process involves a coalescence of distinct crystalline 
clusters 0, HI ■ We find no indication of such behavior in 
our system. At almost all times during the nucleation 
process, the system contains at most one compact crys- 
talline cluster that is of critical size or greater. 

Nucleation thus remains a localized, activated process 
as T decreases toward the HNL, and we do not find evi- 
dence that spinodal-like phenomena influence the nature 
of the liquid state in this range of T. This result is consis- 
tent with a recent study of nucleation in the Ising model 
in which nucleation remained classical even when the nu- 
cleation barrier was less than 8RT [20] . These findings re- 
inforce the strongly kinetic (rather than thermodynamic) 
character of the HNL found here, and the key role played 
by glassy dynamics on approach to this limit. 
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